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Abstract 

We investigate static spherically symmetric solutions of nonlinear mas- 
sive gravities. We first identify, in an ansatz appropriate to the study 
of those solutions, the analog of the decoupling limit (DL) that has 
been used in the Goldstone picture description. We show that the sys- 
tem of equations left over in the DL has regular solutions featuring a 
Vainshtein-like recovery of solutions of General Relativity (GR). Hence, 
the singularities found to arise integrating the full nonlinear system of 
equations are not present in the DL, despite the fact those singularities 
are usually thought to be due to a negative energy mode also seen in this 
limit. Moreover, we show that the scaling conjectured by Vainshtein at 
small radius is only a limiting case in an infinite family of non singu- 
lar solutions each showing a Vainshtein recovery of GR solutions below 
the Vainshtein radius but a different common scaling at small distances. 
This new scaling is shown to be associated with a zero mode of the non- 
linearities left over in the DL. We also show that, in the DL, this scaling 
allows for a recovery of GR solutions even for potentials where the origi- 
nal Vainshtein mechanism is not working. Our results imply either that 
the DL misses some important features of nonlinear massive gravities 
or that important features of the solutions of the full nonlinear theory 
have been overlooked. They could also have interesting outcomes for the 
DGP model and related proposals. 



1 babichev@apc. univ-paris7.fr 
2 deffayet@iap.fr 
3 ziour@apc. univ-paris7.fr 

4 UMR 7164 (CNRS, Universite Paris 7 Denis Diderot, CEA, Observatoire de Paris) 



1 



1 Introduction 



The possibility that gravity is modified at very large, i.e. cosmological, distances 
is motivated not only by the wish to find alternative explanations to the cosmic 
acceleration (or even to replace dark matter [1]), but also because it enables to sep- 
arate what does and does not depend on the dynamics of gravity at cosmological 
scales in the standard cosmological model. Unfortunately, to obtain such modifica- 
tions in a consistent way is very difficult. One simple possibility is to try to give 
a mass to the graviton, and theories of "massive gravity" have attracted recently 
some attention (see [2] for a review). Among those theories, the simplest one are 
what can be called nonlinear Pauli-Fierz theories, i.e. theories made by nonlinearly 
completing the unique consistent quadratic theory for a Lorentz invariant massive 
spin 2, the Pauli-Fierz theory [3J. Such massive gravities share some properties with 
more complicated constructions, such as the Dvali-Gabadadze-Porrati (DGP in the 
following) model [I]. The latter model has also the advantage to produce a late 
time acceleration of the Universe without the need for a non vanishing cosmological 
constant [5] (see however [HI El IB] ) - One of the crucial issues to be answered by such 
models is how to recover metrics sufficiently close to the ones obtained in Einstein's 
General Relativity (GR) to pass standard tests of the latter theory, while at the 
same time having significant deviations from GR at large distances. A first major 
obstacle is related to the so-called van Dam-Veltman-Zakharov (vDVZ) discontinu- 
ity [S], i-e. the fact that the quadratic Pauli-Fierz theory does not have linearized 
GR as a limit when the mass of the graviton is sent to zero. Soon after the discov- 
ery of the vDVZ discontinuity, it was realized that the discontinuity might in fact 
disappear in nonlinear Pauli-Fierz theories [10]. This because a careful examination 
of static spherically symmetric solutions of those theories by A. Vainshtein showed 
that the solutions of the linearized nonlinear Pauli-Fierz theories (i.e. those of sim- 
ple Pauli Fierz theory) were only valid at distances larger than a distance scale, the 
Vainshtein radius Ry, which goes to infinity when the mass of the graviton is sent 
to zero. On the other hand, Vainshtein showed that there exists a well behaved 
(as the mass of the graviton is sent to zero) expansion valid at distances smaller 
than Ry, this expansion being defined as an expansion around the Schwarzschild 
solution of GR. What Vainshtein did not show is the possibility to join together 
those two expansions as expansions of one single non singular underlying solution 
[TT] . In fact, in the context of nonlinear massive gravity, such a joining seems to 
be problematic [T2l [T3] even though it might work for other models where the same 
problem appears, such as the DGP model [TH [15]. Indeed Damour et al. found 
in Ref. [12] by numerical integration of the equations of motion, that singulari- 
ties always appear (at least for the kind of mass term considered there) in static 
spherically symmetric solutions of the kind considered by Vainshtein. A common 
wisdom is that those singularities could be related to the instabilities known to exist 
in nonlinear Pauli-Fierz theories and discovered by Boulware and Deser [11] . Those 
instabilities are also believed [TBI EZ] to be related to the higher derivative operators 
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appearing in a powerful effective description of the scalar sector of massive gravity 
proposed by Arkani-Hamed et al. in analogy with the Stiickelberg "trick" of gauge 
theories [IB] introducing "Goldstone modes" and taking a "Decoupling Limit" (DL) 
in the theory. This limit is aimed at focusing on the strongest self-interactions of 
the model which are seen to arise generically in the scalar sector [18]. It has in par- 
ticular the property to keep the Vainshtein radius fixed and allows to understand 
easily the Vainshtein mechanism as well as the scaling of the dominant terms for 
spherically symmetric solutions in a range of distances between the Schwarzschild 
radius and the Vainshtein radius [HI [161 E] • Our goal in this paper is to investigate 
the success or failure of the Vainshtein mechanism in this DL. It will shed light on 
the possible links between the failure of the Vainshtein mechanism in nonlinear mas- 
sive gravities (as discovered by Damour et al. [12]) on the one hand, the Boulware 
Deser instability [TT] and the Goldstone picture (as introduced in Ref. [18J) on the 
other hand. Clarifying those links matters not only for a better understanding of 
nonlinear massive gravities (theories which are probably of limited interest as far as 
their application to the real world is concerned) but also for the one of more sophis- 
ticated models, such as the DGP gravity, degravitation [Jj5] or the recently proposed 
Galileon [20] ■ For example in DGP gravity, studies of the Vainshtein recovery of GR 
is mostly based on a perturbation theory approach [TH [To] and it matters to know 
if it persists if one takes into account as completely as possible the full nonlinear 
structure of the theory (see e.g. [H |2T]). 

The paper is organized as follows. In the next section (section [2]), we introduce 
with some details nonlinear Pauli-Fierz theories we will be interested in, the Gold- 
stone description of the latter, and the associated Decoupling Limit. This section, 
mostly introductory, will nonetheless also contain some new material concerning the 
Goldstone picture and the DL. We then (section [3]) introduce appropriate ansatze 
for describing static spherically symmetric solutions of the theories introduced in 
the first section, as well as describe with some details the Vainshtein mechanism. 
Next (section H]) we turn to study static spherically symmetric solutions in the DL. 
We first show how to obtain the DL with the ansatz considered, then we study 
solutions of the DL equations of motion. Our main results appear there. We find 
that despite the higher derivative nature of the operator appearing in the DL, there 
exist non singular solutions of the vacuum equations of motion in the DL which have 
the right large distance behaviour. We also show that, depending on the potential, 
this large distance behaviour (at distances larger than the Vainshtein radius) when 
expressed in term of a power serie expansion, might not be enough to select a unique 
solution at small distance. This is in fact a blessing since it leaves more room for 
matching a source at small distance. We also investigate sources and show how our 
non singular vacuum solutions can be extended inside. In addition, we show that 
the small distance behaviour conjectured by Vainshtein is not the only possible one. 
Indeed we find a more generic scaling associated with a zero mode of the nonlinear 
part of the DL equations of motion which also allows the recovery of GR at small 
distance. Interestingly, it allows a working Vainshtein-like mechanism for potentials 



3 



where it was not believed to work. Our investigations (which results are summarized 
in section \5§ also show that the instabilities seen by Damour et al. [12] in the full 
nonlinear case (i.e. not the DL), and hence the failure of the Vainshtein mechanism, 
are not, in contrast to a widespread belief, related to the presence of a ghost (or 
higher derivatives in the DL picture) in the model. 

2 Massive gravity and the Goldstone picture 

2.1 Massive gravity and bigravity theories 

It is well known that the only consistent Lorentz invariant mass term for a spin 
two field Iiab over a Minkowski space-time with canonical metric T]ab takes the 
Pauli-Fierz form given by [3] 

S m = -lm 2 M 2 P Jd^h AB h CD (v AC V BD -V AB V CD ) (1) 

where the kinetic part of the action is obtained by expanding the Einstein-Hilbert 
action at quadratic order in Uab over the flat metric t]ab- I n the equation above m 
is the graviton mass and Mp the reduced Planck mass, given by 

Mp 2 = 8nG N , (2) 

in term of the Newton constant Gn- Wishing to give a nonlinear generalization of 
the Pauli-Fierz theory we have just defined, it seems natural to consider theories 
with two metrics, g^ u and f^ u on a four dimensional space-time^, where one of the 
two metrics, say g fW , will be dynamical, while the other, f^ u , will notd Hence, we 
will consider the action given by 

S = J d 4 xV=g (jfR g + l)j + S int [f, g], (3) 

In the above action, L g denotes a generic matter Lagrangian with a minimal coupling 
to the metric g (and not to the metric /), and Si nt [f,g] is an interaction term with 
non derivative couplings between the two metrics. There is much freedom in the 
choice of this interaction term. For example, the following two possibilities have 
been considered respectively by Boulware and Deser (BD) in Ref. [11] and by 
Arkani-Hamed et al. (AGS) in Ref. [18], 

s$ = -\m 2 M 2 P j 'd* xy /=f H, v H ar (rr-rr) (4) 

S$ = -\m 2 M 2 P j ' d 4 x^-gH, u H aT {g^g^ -g^g"), (5) 

5 For reasons that will appear in the following, we use both capital latin letters and greak letters 
to indicate space-time indices. 

6 Theories of this kind have been considered in the past in the context of strong interactions 
[22] . Note that this is not the only possible way to define a nonlinear completion of Pauli-Fierz 
theory, see e.g. the recent proposal [23] . 
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where f^ u and g^ u denotes respectively the inverse of the metric f^ v and g^, and 
H^ u is defined by 



fiv ■ 



Following the notations of Damour et al. [21], these interaction terms are in the 
form 

SS = -im 2 M| J d'xV^(g, f) = -\m 2 M* J fay/^gV® '(g^f) (6) 

with V^ a \g,f) = y/—g V"( a )(g -1 f) a suitable "potential" density. In this work we 
will consider theories where the potential will not be necessarily one of the 
two above forms (TJHSllll, but will be chosen such that (i) it is a scalar density under 
diffeomorphisms (common to the two metrics) and (ii) that when one expands g to 
second order around the canonical Minkowski metric t\ab as gAB = Vab + ^ab and 
let / assumes the canonical Minkowski form t]abi the potential at quadratic order 
for Hab takes the Pauli-Fierz form ([T]). Note however that many choices can be made 
with those two same properties. The theories considered will then be invariant under 
common diffeomorphisms which transforms the metric as 

9w(x) = df,x' a (x)d u x' T (x)g' aT (x'(x)) , 

Uu(x) = d^{x)d v x' T {x)f aT (x'(x)) , 1 j 

and under which the quantity transforms as a scalar. It is then possible to 
show that it only depends on the matrix g _1 f [23], hence the notation in equation 
(J6]). The equations of motion, derived from action read 

MpG^v = (T, v + I* ) , (8) 

where denotes the Einstein tensor computed with the metric g, T^ v is the energy 
momentum tensor of matter fields, and is the effective energy momentum tensor 
coming from the variation with respect to the metric g of the interaction term Si nt . 
It depends non derivatively on both metrics / and g and is defined as usual as 

A simple, but non trivial, consequence of equations (jSJ) is obtained by taking a g- 
covariant derivative V of both sides of the equations; one gets, using the Bianchi 
identities and the conservation of the matter energy momentum tensor, the con- 
straint 

vz* = (10) 
which the effective energy momentum tensor should obey. 



7 We keep the numbering of the potentials (f4][5]) of Ref. [24] . 
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2.2 The Goldstone picture 



The gauge invariance (J7J) can be used to write the background flat metric / in various 
coordinate systems. Starting from a given gauge, with coordinate X A , and the / 
metric in the form of Jab(X), it might be desirable to change the gauge, but keep 
the change of coordinate explicit in the / metric. Namely, the action is the same as 
action ([3]), but with f^ u (x) now given by the expression 

f^(x) = d IM X A (x)d l/ X B (x)f AB (X(x)), (11) 

while g is kept as g^ v {x). The quantities X A , which then appear explicitly in the 
action of the theory, can be considered as a set of four new dynamical scalar fields, 
which are analogous to the Stiickelberg field used to restore gauge invariance in 
the Proca Lagrangian. Indeed, under a subsequent coordinate change x — > x'(x), 
d^X A {x) transforms as 

d^X A (x)=d^d a X A (x'), 

and hence the quantity fnv(x) transforms as a covariant tensor of rank 2, as it should, 
while Jab is left unchanged by this coordinate transformation. With this in mind 
and the analogy with the Stiickelberg procedure and Goldstone equivalence theorem, 
the initial gauge, where g and / assume the form Qab and Jab is usually called a 
"unitary gauge"[| i.e. one where the Stiickelberg fields X A are gauged away. Note 
that the metric in a non unitary gauge can also be thought as the pullback, 
via the "link field" X A (x), on the space-time manifold 777.4, with coordinates x M , 
of the metric /ab living in an other abstract manifold M.^ with coordinates X A . 
Usually, the unitary gauge is chosen such that in this gauge the metric Jab takes 
the canonical Minkowski form diag(— 1, 1, 1, 1) = T]ab- in the non unitary gauge, 
the action ([3]) is one for a theory with g^ u and X A as dynamical fields. Obviously, 
the equations of motion for g^ v lead to the same equations as in (jSj) where f^ v is 
given in the form pip . Let us briefly turn to the equations of motion for X A . The 
link fields X A only enter the action through the metric f^ u given by equation pip , 
hence the variation of the action with respect to X A , SSi nt (f,g), is given by 

5S int (f l9 ) = -\m 2 Ml J d A x (^l^V^(f,g^j Sfr, 

where 5f^ u is the variation of under X A (x) — > X A (x) + 5X A (x). It is easily seen 
that the latter is given by (see appendix [SJ 

Sf,*, = d^x a Uu + djx a f^ + Sx^dJ^, (12) 

8 In the following, we will always use indices with capital latin letters from the beginning of 
the alphabet A, B,C, ... to denote quantities written in a unitary gauge, like the coordinates X A 
and metrics Jab an d 9ab, and greak letters from the middle of the alphabet, fi, v, ... to designate 
quantities in a non unitary gauge. 
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where 5x^{x) is defined by 



6x»(x) = SX A (x)d A x^(X(x)) 



(13) 



and x^(X) denote the inverse mapping of X J 



[x 



Hence, the variation of / takes 



the form of the Lie derivative Cs x f^u of f along the quantities 5x^ considered as 
a vector field on the space-time manifold. However, the interaction term Si nt (f,g) 
has been constructed to be a scalar under coordinates changes, hence in particular 
under those of the form x M — > x^ + Sx^, with Sx^ defined as above. So one has 



d A x 



5 



Hi/ 



x 



5 



0. 



Using this, and the definition (Q, we get the following expression for the variation 
SS int {f,g) of the interaction term S int {f,g) 



SS int (f,g) 



J d 4 xV=9 (V M Tf ) (d A x°) g av bX J 



(14) 



Hence we see that the equation of motion of X A are equivalent to the Bianchi 
identities ( TTOT) provided that the mapping X A (x) is invertible. Note that above, we 
have lifted and lowered indices with g^ v . This will be the case in all the following 
except for what concerns the metric f^, the indices of which are moved with / 
itself. Note also that in the above derivation we have not used the fact that Jab is 
the canonical Minkowski metric tjab and hence, it also holds when it is not the case. 



2.3 Strong coupling and decoupling limit in the Goldstone 
picture 

The authors of Ref. [18] have further developed the above mentioned analogy be- 
tween X A and Stiickelberg fields. Indeed, following [18], one can do a "Goldstone 
boson" expansion of the action ([3]) around a unitary gauge. Indeed, considering 
some background solution for g^ v (defined as g® u ) and X A (x) defined (the metric 
f ab being kept fixed) as 

X A (x) = 6 A x», 
Ref. [18J introduces the "pion" fields n A as 

X A (x) =X A (x) + 7t A (x), (15) 
and further does a " scalar- vector" decomposition of the n A in the form 

n A (x) = f AB (A B + d B 4>). (16) 

Note that this notation is a bit problematic if one is interested in the structure 
of the theory beyond linearized equations of motion. Indeed, first neither the n A 
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defined as above, nor and A B are tensors on the Manifold M.^ (and one does 
not see the need for the metric f AB above). Second, given the definition (j!5p . it is 
natural to assume that is a function of x^ (not of the X A ) and hence one does 
not see why it is differentiated above with respect to X B . In fact, if one takes the 
definition (fT6l) literally and write 8b<P = dsx^d^, one sees that dsx^ is expressible 
formally as a serie of n A (using ffT5]) ) and hence also of (f>. This will generate terms 
at nonlinear order which will be different from the ones one gets using the naive 
expression <9s0 = S^d^. The same will be true if one considers metrics Jab with 
non trivial dependence in the coordinates X A (so that this later problem does not 
arise when f^ B is taken to be the canonical flat metric, but will e.g. if one chooses to 
parametrize Minkowski space-time in a non trivial way). To avoid those difficulties, 
we will write, instead of ( 1T61) (and in fact this seems to be what is done implicitly 
in Ref. [18] to deal with nonlinear order) 

ir A (x)=6 A (A< l (x)+r)^d u( f>). (17) 

If one inserts the decomposition (1171) into action (jHJ), and expands around flat space- 
time writing g^ u = r]^ u + h^, we obtain an action for the dynamical fields h^x), 
A^{x) and 4>(x). Since A^(x) and 4>(x) only enter in the metric f^u, via expression 
pip , the only term in action ([3]) which depends on A^(x) and <f>(x) is the interaction 
term Si nt [f, g] , and one has from ffTTl) (where no term has been neglected in the 
expression below) 

= h^u - d^Ay - dyA^ - 2<9 M <9 i ,0 

dp L A a d u A <T - d„d°(j) (18) 

-d v A a dpd a 4> - d„A a d v d a (j). 

Inserting this expression into S[f,g], and keeping the lowest order in h^(x), A fl (x) 
and <f>(x), one obtains the following first non trivial terms (upon integration by part) 

M 2 f ( 

S = -^ d 4 x{2h^d^d u h - 2h^d v d a h° + h^Dh^ - hUh 

+ m 2 [h 2 - h^ v - F^F^- A(hdA- h^A»)-4(hn<J) - h^d^cj))] } 

+ 

where h = h flu r] tJjU , dA = d^A^ 1 , T^ u is the matter stress-energy tensor, and indices 
are moved up and down with the metric T]^. The peculiarity of the above expression 
is that while A 11 acquires a standard kinetic term, <fr does only get one via a mixing 
with h^ u [18] . this being entirely due to the structure of the Pauli-Fierz mass term 
([TJ). A may to demix and is to do the shift [18] 

hfa, = hfiu ~ fn 2 !]^. 
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The quadratic interaction action then becomes 



M 2 m 2 f r~ - 
S^t = — |— / d 4 x{/i 2 - - F^ U F^ V - 4(hdA - V,3M") 

+ 6m 2 <f>(D + 2m 2 )(f) - h<p + 2<f>dA J. 

The interactions between 0, A and /i^ can be canceled by adding an appropriate 
gauge fixing to the action (see [25J ) . Following again [18], we can obtain canonically 
normalized fields 0, A and h^ u by defining 

h^v = Mph/j, u , 

A" = MpmA^ (19) 
= M P m 2 (j). 

Doing so, and expanding the action in 0, A and h^ v one sees using ( TTSil that has 
in general cubic self interactions suppressed by the energy scale 

A = (m 4 M P ) 1/5 . (20) 

When those interactions are present, they are the strongest interactions among the 
fields 0, A and h^ u in the limit where m <C Mp, besides quadratic, cubic and quartic 
non derivative interactional Indeed, using the expansion fl 1 8H . it is easy to see that 
those interactions, coming from the interaction term ([6]), scale as 

A^SSr 8 *"-^* 1 (9A) k2 (dd$) ^ ^ (21) 
where ki, k 2 , k 3 and fc 4 are integers and the expressions A fcl fc2jfc3ifc4 are given by 

4fc 1 +3fe2+2fc3+4fc 4 -6 
/ Mp \ 5(fc 1 +2fc 2 +3fc 3 + fc 4 -4) 

^kxMMM = ^ \ ■ ( 22 ) 

\ m J 

Whenever ki+2k 2 + 3k 3 + fc 4 = 4 (which can only happen for ki+k 4 = 4, k 2 = 0, k 3 = 
or for ki+k 4 = 2, k 2 = 1 and A; 3 = 0), the coefficient in front of h kl (dA) k2 (<9<90) fc3 fe4 
in Eq. (|2ip is equal to the dimensionless expressions m 2 /Mp or m/Mp and the 
corresponding interactions include in particular quartic non derivative couplings 
between and h. In all other cases, Ak ly k 2 ,k 3 M has the dimension of an energy. 
It is equal to m for quadratic interactions and m? /Mp for cubic non derivative 
interactions. Hence, besides quadratic, cubic and quartic non derivative interactions, 



9 An appropriate choice of the interaction term Si n f[f, g] can remove cubic (and others) self 
interactions of <j> [18] . 

10 We include here and in the following with those interactions, a peculiar set of cubic interactions, 
with a derivative coupling to A, corresponding to k\ + fc 4 — 2, hi — 1 and k% = (see below), 
which will play a similar role. 
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it is easily seen that the strongest interaction among the canonically normalized 
fields is the cubic derivative self interaction of <fi suppressed by the scale A = A ,o,3,o- 
All the other interactions are suppressed by a scale larger or equal to A 0j o,4,o> with 
Ao.0,4,0 — Mp 4 m 3 ^ 4 [T8| [T6l [25]. In other words the exponent in the right hand side 
of the above formula (j2~2~i) is bounded below by 1/20 as soon as {^1,^2,^3,^4} 7^ 
{0, 0, 3, 0} (not considering the above mentionned non derivative quadratic, cubic 
and quartic interactions). Hence, for m -C Mp the scale A4 is stricly larger that A, 
while the non derivative cubic and quartic interactions are much smaller than the 
quadratic mass terms (for small expectation values of the fields). This indicates that 
there is a regime where the theory considered, and hence also its solutions, is well 



approximated_J by retaining only the quadratic action and the cubic self interaction 
of 0, as noted in particular in Ref. [TB]. This regime can be extended to arbitrarily 
high energy scale 12 ! (or arbitrarily small distances) by choosing a sufficiently large 
ratio Mp/m, as can be seen from the relation (122]) . In fact, one can take a decoupling 
limit, that suppresses all the interactions but the cubic derivative self interaction. 
This limit is defined as 

M P -> 00, 
m ~ * 0, 

( 23 ) 

A ~ constant, 
T^/Mp ~ constant. 

In this limit, all the quantities A^^MM (with {ki, fc 2 , k 3 , k±} 7^ {0,0, 3, 0}) are sent 
to zero or infinity. The action one is left with for is of the form! 13 ! 



2 



3 7 



1 



d*x { -0D0+ — 



a (D0) 3 + (3 (□<; 



M, 



-T<b 



(24) 



where a and (3 are numerical coefficients that can be adjusted at will by choosing an 
appropriate interaction term S in f[f, qfn. and T is the trace of the energy momentum 
tensor. For example, the BD potential (j4]) leads to a = —f3 = —1/2, while the AGS 
potential (jSJ) leads to the opposite case a = —(3 = 1/2. The equation of motion 
deriving from this action is 



3n □ \n4>) 2 + (3D 



+ 2(3d,d u □ 



-T . (25) 



This equation of motion is fourth order which signals generically ghost propagating 
degrees of freedom. In fact, one can argue [TB1[T7] that one can see this way a generic 



11 Of course the discussion here is a bit loose, it will be made more precise when dealing with 
spherically symmetric solutions. 

12 We do not consider here, as done e.g. in [TS], the issue of the quantum corrections to the 
theory and only discuss it from a purely classical perspective. 

13 While h A become free. 

14 Note that in general, the cubic term for cf> is given by some linear combination of the three 
terms (D0) 3 , □</> <p,^i> 4>'^ 1 ' and ,MQ but an integration by part can always be used to 
reduce the number of independent terms to two, as shown in Eq. (|24|) . 
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property of massive gravity once discovered by Boulware and Deser [UJ: namely the 
fact that nonlinear massive gravity propagates at nonlinear level one more degree of 
freedom than linear Pauli-Fierz theory, the energy of the extra propagating mode 
being unbounded below. This is a major obstacle to the possibility to consider 
nonlinear Pauli-Fierz theories defined as above as a realistic theories (see however 
|26j), but, again, our aim is here to use this theory like a toy model to study the 
Vainshtein mechanism, and not to advocate for a realistic use of it. 

3 Static Spherically symmetric solutions 

In this section we introduce our framework to look for static spherically symmetric 
solutions of massive gravity. We also describe the Vainshtein mechanism. First, in 
Sec. 13.11 we present ansatze for the metrics and discuss possible coordinate choices. 
In what follows we will mainly be interested in the bi-diagonal ansatze, i.e. those 
where both metrics can be put simultaneously in a diagonal form. It turns out that, 
in this convenient coordinate choice is so-called the A, //, v gauge (following 

the terminology of reference [12]) that will be introduced. The equations of motion 
in this gauge (together with the Bianchi identity) form a system of ordinary differ- 
ential equations for the three functions A, v and fi to be determined. This system is 
strongly nonlinear and cannot be integrated analytically. However, different asymp- 
totic regimes can be studied separately. First, one expansion can be made in the 
Newton constant, Gn- This expansion will be found to be valid far from the source 
(i.e. for distances R 3> Rv, where Ry is the Vainshtein radius that we will intro- 
duce). At lowest order, it does not match the similar expansion that one can make 
in General Relativity, this being reexpressed in the form of the vDVZ discontinuity. 
Vainshtein conjectured that close to the source the General Relativity solution can 
be restored via the effect of the nonlinear corrections. This is the essence of the 
Vainshtein mechanism explained in section 13.21 In the Vainshtein original proposal 
[TO] , a scaling with distance is proposed for the first correction to the Schwarzschild 
solution close to the source (i.e. for R <C Rv)- Here we show that there is another 
possible scaling that will play an important role, as will be discussed in section HJ 

3.1 Ansatze and coordinate choices 

The most general static spherically symmetric ansatz for solutions of the theory 
defined by action ([3]) takes the form 



cj ^^jdoc doc 
f fiiydcc doc 



V 



J{r)dt 2 + K{r)dr 2 + L(r)r 2 dtt 2 , 

C(r)dt 2 + 2D(r)dtdr + A(r)dr 2 + B(r)dtt 2 , 



where dfl 2 



is the canonical metric of a unit 2-sphere 



dtt 2 = d9 2 + sin 2 6dtp 2 . 
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The above ansatz can further be simplified by setting L be one, by a suitable choice 
of the radial coordinate r. However, in general, it is not possible to put both metric 
at the same time in a diagonal form by a coordinate change without imposing restric- 
tions on the solutions, this is because we only have one diffeomorphism invariance 
(J7J) at hand. Different cases have been considered in the literature and some solu- 
tions are explicitly known in the case when metrics are not both diagonal [27] [12] . 
In this work we will only consider bi-diagonal cases where the non dynamical metric 
is parametrizing a Minkowski (non dynamical) background space-time. It will turn 
out useful to use different type of gauge and ansatze. A first gauge choice (called 
the "a,b,c gauge" in [12]) is defined by metrics in the following form 

g AB dx A dx B = -J(r)dt 2 + K(r)dr 2 + L(r)rW 

f AB dx A dx B = -dt 2 + dr 2 + r 2 dQ 2 ' A> ' 

A simple coordinate change X A (X') 

z = r cos 9, 

x = r sin 9 cos ip, 

y = r sin 9 sin (p, 

puts of course Jab m the canonical Minkowski form t]ab (we will not need the 
corresponding expression for g) 

f AB dx A dx B = -dt 2 + dx 2 + dy 2 + dz 2 . (27) 

Note that both gauges f)26p and (j2"7|) have in common that all the unknown functions 
(J, K and L) are put in the metric g and hence, in agreement with our definitions 
of section 12.21 we can consider those gauges as unitary. A non unitary gauge, which 
has some advantages, is the one where one of the unknown function J, K, L is put 
into the expression of the non dynamical metric, given by 

g^dafdx" = -e v{R) dt 2 + e x{R) dR 2 + R 2 dQ 2 , 

f^dx»dx v = -dt 2 2 e-^dR 2 + e-^R 2 dQ 2 , (28) 

where here and in the following a prime denotes a derivation with respect to R. The 
relation between gauges (|2T)j) and (12"%1) is given by 

ife-M(*)/2 = r? ( 29 ) 

while the relation between the functions appearing in the metric coefficients are 

J(r) = e u{R) 



K(r) = e^ H >e 
L(r) = e» {R) 



2 

a; /.'),/-(/.') 
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Following reference [12] . we will call the gauge (1251) the gauge. It has the 

advantage that the g metric can readily be compared to the usual Schwarzschild 
metric of standard General Relativity. The coordinate change defined by 

Z = R cos 9, 

X = R sin 8 cosip, 

Y = R sin 8 simp, 

puts then the metric g and / in the form 

g^dx^dx" = -e u{R) dt 2 + dX 2 + dY 2 + dZ 2 



(30) 



x r - 

R 


= rd x R, 


r 

Y R 


= rdyR, 


z r - 

R 


= rdzR- 



+ ( — — 2 — J {xdx + YdY + zdzy 

f^dafdx" = -dt 2 + e~^ R \dX 2 + dY 2 + dZ 2 ) 

+ (^-VL + Vpj e -m (XdX + YdY + ZdXf 

with R 2 = X 2 + Y 2 + Z 2 . We obtain easily the relation between coordinates = 
{t, x, y, z} of the unitary gauge ( 1271) and coordinates {x^} = {t, X, Y, Z} of the non 
unitary gauge ( 1301) as (the time coordinate t being the same in the two gauges). 

x= Xe~^' 2 

y = Ye-^ R) ' 2 = Y^ = rd Y R, (31) 
z = Ze-^l 2 

R 

The Jacobian J of the transformation (|3ip . I^X" 4 !, is given by 

J = e-l^) (l - . (32) 

It does not vanish except possibly on a sphere of radius R where 

= 1- (33) 

Obviously, a generic function \i would not define via Eq. ( 13T1) a one to one mapping 
over the whole Minkowski space. To obtain such a mapping, a sufficient condition is 
that r(R) is a strictly monotonous function of R which maps the real positive half 
line to itself with an everywhere non vanishing Jacobian. In this paper, because we 
will deal with the Decoupling Limit, we will not be careful about what is happening 
at the origin r = 0. However, a necessary condition on fi is that r(R) maps the 
origin to itself, and hence \x should verify 

lim Re~^ R)/2 = 0. (34) 
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Note that the monotonicity of r(R) is automatic, provided that (|33|) never holds. 
Given a particular solution for /i, A, and v, the above conditions could prevent the 
solution to be interpreted as a correct solution of the equation of motion of the 
theory (j3J). 

3.2 The Vainshtein mechanism 

3.2.1 A short introduction to the vDVZ discontinuity and the Vainshtein 
mechanism 

The (quadratic) Pauli-Fierz theory, with a mass term given as in Eq. (JTJ), is known 
to suffer from the van Dam-Veltman-Zakharov (vDVZ) discontinuity, i.e. the fact 
that when one lets the mass m of the graviton vanish, one does not recover predic- 
tions of General Relativity. E.g., if one adjusts the parameters (namely the Planck 
scale) such that the Newton constant agrees with the one measured by some type 
of Cavendish experiment, then the light bending as predicted by Pauli-Fierz theory 
(and for a vanishingly small graviton mass) will be 3/4 of the one obtained by lin- 
earizing GR One way to see this is to consider solutions of equations of motion 
(jHJ) which are static and spherically symmetric and which would describe the metric 
around a spherically symmetric body such as the Sun. To do so, using the ansatz 
(|28|) is especially convenient because in this form the g metric can be easily compared 
with the standard Schwarzschild solution. If one tries to find a solution expanding 
in the Newton constant, as we recall in subsection 13.2.31 one finds immediately the 
vDVZ discontinuity appearing in the form of a different (m independent) absolute 
value of the coefficients in front of the first non trivial correction to flat space-time 
in g tt and g^R components (neglecting the Yukawa decay by assuming the Compton 
wavelength of the graviton is much larger than other distances of interest). How- 
ever, as first noticed by Vainshtein [10J, the computation of the next order correction 
shows that the first order approximation ceases to be valid at distances to the source 
smaller than a composite scale, the Vainshtein radius defined by 

R v = (m- 4 i? 5 ) 1/5 , (35) 

where Rs is the Schwarzschild radius of the source. This Vainshtein radius obviously 
diverges when one lets m go to zero and in fact is much larger than the solar system 
size for a massive graviton with a Compton wavelength of the order of the Hubble 
radius. Hence, one can not conclude that massive gravity is ruled out based on solar 

15 The fact it is smaller is easy to understand: the essential difference between Pauli-Fierz theory 
and linearized GR comes from an extra propagating scalar mode present in the massive theory. 
This mode exerts an extra attraction in the massive case compared to the massless case. Hence, 
if one wants measurements of the force exerted between non relativistic masses to agree, the 
coupling constant of the massive theory should be smaller than that of the massless theory. But 
light bending is blind to the scalar sector - because the light energy momentum tensor is traceless. 
Hence, provided the two theories agree on the force between non relativistic probes, the massive 
theory would predict a smaller light bending than the massless one. 
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system observations and results of the original works on the vDVZ discontinuity [9] . 
Vainshtein also showed that an expansion defined around the standard Schwarzschild 
solution can be obtained (as recalled in subsection 13.2.41) that is well behaved when 
the mass of the graviton is sent to zero, opening the possibility of a recovery of GR 
solution at small distances R of the source. Indeed, the domain of validity of this 
second expansion was shown to be R <C Ry. Moreover, the correction found by 
Vainshtein to the Schwarzschild solution are non analytic in the Newton constant 
which could have explained the failure of the attempt to obtain a solution expanding 
in the Newton constant. This is the aim of this section l3~2l to give more details on the 
so-called Vainshtein mechanism, i.e. the possibility of a non-perturbative recovery 
of solutions of GR, which is the bulk of the studies of this work. 

3.2.2 Equations of motion in the A,/i, v gauge 

Let us first consider the gauge (128|) . In this gauge, the equations of motion (JE]) read 



where the source energy momentum tensor T " is assumed to have the perfect fluid 
form 





8nG N (T 9 RR + Pe x ) , 



(36) 



V 



diag(-p, P,P,P), 



with total mass 



M 




The matter conservation equation reads 



P' 




(37) 



while the Bianchi identities have the only non trivial component 



1 1 



(38) 



m 2 Mp R 



0. 
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Expanding equations (1361) and (1551) in power of A, v and p, we find 



^ + ^ + ^%,A) + Glf W) = -^!(A + 3// + V) + 87rG' i vp 

+m 2 Q u (jy, A, /x) + vMG N p 
+m 2 C? t (v,\,p) 

+vMG N p(e» -1-v), (39) 

+m 2 Q RR (u, A, p) + X8ttG n P 
+™ 2 C^?(^A,/i) 

+v8irG N P(e x - 1 -A), (40) 
^-^-^(^A,/i)-C>(z/,A,/i) = 0, (41) 

where Cr# , G^ RR and represent respectively the quadratic (in power of A,/i 
and u) part of G/y? and the Bianchi identity; G tt , G RR and represent 
respectively the cubic and higher part of G u , Grr and the Bianchi identity. The 
expressions of some of those quantities can be found in appendix [B] 

3.2.3 Expansion in the Newton constant 

We can first look for solution of the system (I39H41I) expanding the solution into 
powers of the usual Schwarzschild radius Rs of the source (or into the Newton 
constant). More explicitly, we expand as in 

A = A + Ai + ... 

v = V + l>i + ... 

p = po + pi + ... 

where \i,Ui,pi are expected to be proportional to G 1 ^ 1 , and we consider a regime 
where Aj + i <C Aj, z/ i+ i <C Vi and pi + \ <C p^. 

At linear order, far from the source, the equations of motion reduce to 



R R 2 2 



+ -7^ = — — (Ao + 3/i + Rp'o), 



5-F = T^ + ^- (42) 

Ap _ v'q 
R 2 ~ 2R' 

Notice that the expansion of the Bianchi identity does not contain a term linear in 
p, as can be seen from equation fj4Tl) . this being due to the peculiar structure of 
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the Pauli-Fierz mass term. The exact solution of the system (|42|) can be found in 
Ref. [12]. However, we are interested here in the limit where Rm <C 1, and want to 
simplify the linear system accordingly. One can see from the third equation that Ao 
and z/o are of the same order: Ao ~ vo, while the first or the second equation indicates 
that /io ~ Ao/(mi?) 2 Ao, v§. As a consequence, the system can be simplified to 

A'o . A m 2 , 

-ft+tf = -^-( 3 /^o + ^o)> 

Vq A 2 

= m/i °' 
Ap_ Vq_ 
R 2 ~ 2R 7 

with the scaling (valid for mR <C 1) 

/io > A ~ v Q . (43) 
One can extract from the above system an equation for /io reading 

3/i + Rfi'o = 0, (44) 

which can easily be solved to get 

\ _ Ci I Ci 

Ao " 2R' ~ to ~ T^kJ 2 2R • (45) 

where C\ is a constant of integration which is expected to be proportional to 
and has to be fixed by matching to the source. 
At second order, we solve 

^ + Jj2 + G t?\"o,\o) = -^(Ai + 3/xx+iVi) WQ«KAo,//o), (46) 

- ^ + G^\u , \ ) = — + +m 2 g rr (z/ , A ,/x ), (47) 

it 2 = ^ + £ &KAo,/*o). ( 48 ) 

Assuming that v\ and Ai are of the same order (in line with u ~ Ao), and using the 
scaling (j4"3~|) we find from fT4H|) the scaling 



\ 2 Rs ^ 1 Rs 

1 1 Po i? (mR) 4 R 

As a consequence the above system (I46H48I) reduces to (in the mR <C 1 limit) 

Ai Ai Tfi 2 ... .. * s , 

= ™V (50, 

F - + < 51 > 
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with Q{[1q) = Qb(0,0, fi ). For the interaction terms of Eq. (j4]) and ([5]), we find 
respectively (see the expressions given in appendix IBI) 

g«0.) = ^ + ^ + (52) 

QW(A) = -Q ,2, W. (53) 
In the general case one can show (cf. appendix [B]) that Q(fi) is given by 



{ 5 3 

+[3 I 10//// + 5iV 2 + yRw" + 2 R2 ^'^" 



(54) 



where a and /3 are numerical constants depending on the interaction term Si nt [f, g\. 

We can eliminate Ai and i>\ from the above system ( 1491I5T1) to get an equation 
for fix 

3 

3Q(yWo) + rQ'(fio) = ~^nfi 2 (3/ii + ityi'i) • 
Hence, gathering the linear and quadratic terms in the Gn expansion, we find 

* - ^4w +o(j|1 (56) 

where m, mi and Zi are order one dimensionless quantities that depend on the 
details of the quadratic term Q(fio), and Rs is the Schwarzschild radius defined 
as usual from the "Cavendish" Newton constant Gn (defined as the one which is 
measured in some Cavendish-like experiment, assuming the Newtonian potential is 
given at leading order by the first term in the right hand side of equation fl55l) ). To 
obtain the expressions (15511571) above, we have fixed the integration constant C\ such 
that the Newtonian force between non relativistic pointlike bodies matches the one 
obtained in Newtonian theory (assuming those bodies are separated by a distance 
much larger than the Vainshtein radius and much smaller than the graviton Compton 
wavelength). As we said above, this requires defining the Newton constant Gn as 
4/3 of Gnl3 as defined by equation Eq. (j2j). We also see the vDVZ discontinuity 
in the fact that the coefficient in front of Rs /R in the first term on the right hand 
side of (156j) is not equal to one (as it would be in GR). Finally, we see that this 
expansion is only valid for R ^> Ry where Ry is the Vainshtein radius defined as in 



16 So that we have also Rs = 4/3i?s, where Rs is the Schwarzschild radius used in the rest of 
this paper. 



18 



3.2.4 Small R expansion 

At small R, R <C Ry, following Vainshtein's idea, one is looking for an expansion 
around solutions of usual General Relativity. In other words, one looks for an expan- 
sion in power of the mass of the graviton (squared) m 2 and expands the functions 
A, v, \i as in 

oo 

f(R) = Y,™ 2n fn(R), 

n=0 

where the f n are m-independent coefficients. By definition, the lowest order expres- 
sions of A and v are given by the Schwarzschild form 

^--- h ('-fi4(i) V - (58) 

Keeping the lowest order term in R$/R, Xq and v$ are simply obtained by solving 
the vacuum linearized Einstein equations 

^> + ^ = 
R R 2 



v'o A 



0. 



R R 2 

The function /x is found using the Bianchi equation 

£ = Is + «<"•>• 

where only the lowest order terms have been kept. Inserting in the above equation 
the GR solution (|58p . we find that one should have 

« = 7^§, (59) 

where the lowest order in Rs/R has been kept. Assuming, following Vainshtein, 
that [Mq can be expanded as a power of R, we find that the first non trivial term 
should be of the form 

fio = M ^R S /R } 

> A ,v 0) (60) 

where Mq is a pure number. Hence, there is a first potential obstruction to the 
success of the Vainshtein mechanism [12], namely the left hand side of equation 
(159]) must be positive definite, which is only possible for particular quadratic terms 
Q(Ho). E.g. the interaction term (j3j) leads via equation fl59l) to an imaginary M 
and hence, for this potential the scaling proposed by Vainshtein does not work (see 



19 



however below). In contrast, from the interaction term (JSJ) (and hence equation 
(I53"|) ). one finds the real value M = ± a/8/9. 

The next order for Ai and /ii is given by solving equations ( 13911401) . In the limit 
where 

Rs < R < m -1 , 

taking into account the scaling (|60|) . those equations simply reduce to 

A'i Ai m 2 , „ , . 

i + ^ = -^-( 3 ^o + M) 
Ai 



"77 77J = m 2 /i , 



while /ii is obtained solving the Bianchi identity fT4"Tj) which reduces, in the limit 
considered, to 

where in the right hand side of the above equation only the term linear in /io and 
/ii are kept (cross products). Eventually, we obtain an expansion of the form found 
by Vainshtein given by 



v 



~ + N, (mR) 2 sJ^ + O (m 4 ) , (61) 



A = ^ + L x (mR) 2 \^- + (m 4 ) , (62) 




/i = M \ -± + M 1 (mRy + (D(m A ) , (63) 



where again Mo, N±, L\ and M\ are order one dimensionless numbers. This expan- 
sion makes sense only if m 2 ^R 2 ^^J [1 + (^f )] < [l + O (^f )] , i.e., if 

i* < = (m-*R s ) 1/5 . 

Note however that there is another possible expansion at small R that has not 
been discussed previously in the literature and for which one also recovers the GR 
solution. Indeed, the operator Q, appearing on the left hand side of equation fl59l) . 
can have a zero mode in the form of a power law, fiQ oc R p , and one can find an 
expansion for starting with this zero mode. E.g. in the cases of interaction terms 
(|4|) and ([5]), where the corresponding operators Q are given by equations ( 152|) and 
( 153]) . one finds p = —2, and hence in those cases, as we will discuss in more details in 
the next section (together with a discussion of the most general case) it is possible 
to find a solution to equation (1591) in the form 

= (mRsf 5 ( A 2 - A log(R/R v )\ + O (JL\ , (64) 
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where s = — 1 for the potential (jl]) and s = +1 for the potential (jSJ), and A is a 
pure number. Such an expansion will in fact, as will be shown in the next section, 
turn out to be the most general one found by numerical integration of the equations 
of motion in vacuum. 

3.2.5 Matching small R and large R behaviours 

It was noticed immediately after the seminal work of Vainshtein that there was 
no warranty that one could match the solution found by Vainshtein in the small 
R limit, given by equations ()61H63p . to the one obtained in the large R limit, and 
given by equations (15511571) [TT] . More recently, numerical integration of the full 
(nonlinear) system of equations (I39H41P has shown indeed [12] (see also p3]), that 
integrating inwards from the large R behaviour (I55H57I) or outwards, from the small R 
behaviour, always results in singularities appearing at finite R (provided one insists 
upon the solution to be asymptotically flat). In fact the exact reason for which those 
singularities arise has not been clarified so far in the literature and it is one of the 
main purpose of this work to reexamine this question in the light of the Goldstone 
formalism in the decoupling limit of Ref. [18]. The general understanding is that 
one can link those singularities, and hence the failure of the Vainshtein mechanism, 
to the ghost-like nature of the extra degree of freedom discovered by Boulware and 
Deser [TT] to be present in the nonlinear theory, and, as recalled above, appearing 
in the Goldstone formalism in the form of higher derivative operators [16j [TT]. 

4 The decoupling limit in spherically symmetric 
solutions 

In this section we study static spherically symmetric solutions in the decoupling 
limit (DL in the following). We first show (section 14. 1 . lj) how to obtain a DL in 
the case of the static spherically symmetric ansatz (|28p . We then show (14.1.21) that 
this DL is in fact analogous to the one introduced in section 12.31 and we show the 
scalar field is closely related to the metric function /i of the gauge (|28|) . We 
give the explicit link between the two and show that in the decoupling limit, both 
fields obey the same equation of motion. The DL allows to decouple the scalar 
degree of freedom, which appears due to the presence of the potential term, from 
the degrees of freedom also present in GR. In this limit, the study of solutions is 
greatly simplified, in particular, instead of three equations of motion (I39ti4ip there 
will be only one (nonlinear) equation for the "gauge" function fi (or equivalently for 
(f)). This is consistent with the fact that by extracting the DL the GR degrees of 
freedom become decoupled from the scalar degrees of freedom. 

On the other hand, one can hope that the DL will contain the main features 
of the full system, since it isolates the "troublesome" scalar degrees of freedom. 
For example, the large-i? asymptotic solution and the Vainshtein solution of the 
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full system have DL analogues. On the other hand, as was recalled at the end of 
the previous section, the Vainshtein mechanism was found to fail for the theory 
considered here, by integrating the full nonlinear system of equations ()39H4ip . Thus 
a natural question to ask is whether one can already see this failure in the DL. In 
the section 14.1.31 we discuss briefly the expected range of validity of the DL. 

The main results of this section are obtained in section 14.21 where we solve nu- 
merically the DL equations of motion. After a discussion of boundary conditions 
at infinity (subsection 14.2.11) . we first introduce a simple linear differential equation 
which shares some crucial properties with the one obtained in the DL (subsection 
I4.2.2p . Then we discuss solutions outside the source and inside the source in the 
simplest cases of potential (j4]) and ([5]) (subsections 14.2.31 and I4.2."4|) . Last, we turn 
to discuss the most general case (subsection (14.2.51) . 

4.1 Extracting the limit 

4.1.1 Rescaling and limiting equations of motion 

First, let us note, as should be clear from the derivations presented in section 13.21 
that both the expansion in the Newton constant (in the limit where R <C m _1 ) 
as given by equations (I55H57I) and the expansion in the mass of the graviton (in 
the limit where R 3> R s ) as given by equations (I61ti63|) (as well as (164")) ) are left 
unchanged if one first does the rescaling on the metric functions, 



and then takes the DL as it is defined in fl23|) . More precisely, plugging the rescaled 
functions z>, A and fx into the equations of motion (13911411) . and assuming those func- 
tions to be bounded, one sees that in the DL (T2"3"l) . the nonlinear system of equations 
(ESmD collapses to the much simpler one 



v = M P v, 
A = M P A, 
jl = m 2 MpLi, 



(65) 



and on the components of the energy-momentum tensor 



P = p/Mp, 
P = P/Mp, 



R + W 
R R 2 



-(3fx + Rfl') +p 



R 2 




(66) 
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where A is denned as in ( 120]) . So the only non linearities which remain in this 
limit are those corresponding to the piece of the Bianchi identity quadratic in p. 
Note that we have set the pressure P to zero in the above equations (IBTij) . since this 
is a direct consequence of taking the decoupling limit in the matter conservation 
equation (J3"7j) : 

v' 

p> = -— — (p + p) _> when M P -> . (67) 
2Mp 

It is easy to extract from the above system fl66l) a single equation obeyed by p 
reading 

i [6Q(p) + 2J2g(/2)'] + ~£ + ~i? p' = p . (68) 

As we will see in the next subsection, one can map this equation to the equation of 
motion found for <ft, equation and hence there is a very clear relation between 
the decoupling limit as we just defined, obtained from the equations of motions in 
the A, p, v gauge, and the one discussed in the original reference [18j as summarized 
in section 12.31 Before turning to discuss in detail this relation, let us first note that 
equation (1681) can easily be integrated once leading to the first integral 

^W) + \* = § + ^ J*** ~p{ k ) R2 (69) 

where is an integration constant. Notice that outside the source, the integral on 
the right hand side of equation fl6"9"|) gives a constant and hence this simply results 
in a shift of the integration constant C 2 to yield 

|a W) = (™) 

As usual, the integration constant C should be fixed by matching to the source. 
The two limiting regimes discussed in section 13.21 are easily recovered from this 
equation. The lowest order term in the expansion into the Newton constant (valid 
at R Ry) corresponds to keeping only the linear order in p in the left hand side 
of (17D1 . recovering equation (PHI) . The lowest order term in the expansion in m? 
(valid at R Ry) corresponds to keeping only the quadratic term Q(p) in the left 

hand side of ATP]) , recovering the scaling p ~ O yy/Rg/Rj found in equation (!63j) . 
Assuming the validity of the Vainshtein recovery mechanism, one can easily fix the 
(shifted) integration constant Co for a point like source, by using the known form of 
the functions v and A at lowest order in the m 2 expansion which is given by their 
expressions in standard General Relativity (at lowest order in Rg/R). This fixes 
that, in this limit (see e.g. equation (|59|) ). 

-|i? 3 Q(p) - M P Rg, (71) 

and hence, 

Co = M P Rg . (72) 
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Note that the above equations (I69H70P take their most simple expressions when Q 
is given as in equations (I52H53I) . i.e. for the interaction terms (jlj) and ([5]), in which 
case e.g. equation (1701) reads 



A5 + 



I + -p = — 



R J 



R3 



(73) 

where again s = — 1 for the BD interaction term (T4|) and s = +1 for the AGS mass 
term ([5]). For future use, and for numerical integration, it turns out convenient to 
introduce the dimensionless quantities 





= Rj Ry , 


Pa 


A 

= Air — -p , 


HO 


= cT 2 /i , 


v(0 


-4 

= a v , 


«(0 


= a -4 A , 



where we used the dimensionless parameter 

a = R v m = (R s m) 1/5 . 

The functions u and v are just the dimensionless functions A and v associated with 
the g^ u metric, while w corresponds to /x, and £ is the dimensionless distance R 
expressed in unit of the Vainshtein radius Ry. In term of these quantities, the 
system ( Itjtj]) reads 



u u 

I + e 

V u 
u 

e 



-~(3W + £W) + p a 



w 



where a dot denotes a derivative with respect to £. Eliminating u and v from the 
above system, we get the analogous of Eq. ( 1701) for the rescaled variable w 



2 Q(w) + \w = p, 



(74) 



where Co is defined by 



c 



Co 



Fixing the integration constant as in (172]) . leads to c = 1, and the equation for w 
now reads 

2 QH + lw = i (75) 



24 



which will be assumed from here-on and until the end of this article. The general 
form of Q can be read from equation (1541) . R being replaced by £ and \i by w. It is 
given by 

1 I „ ( ' i . .. 3 .. _ . 2 Gww\ 
-WW H — WW + 2w H ; — 



Q(w) = — < 3a I -ot H — ww + 2w + ( , 
2 I \2 2 s / 



(76) 

* . .. 5 .. r . 2 KtemA 

+ p I yOT + -TO + OW H I 



E.g., in the simplest case of the potentials leading to equation (J73J), one has a+(3 = 
and equation (175]) reads 

We now turn to compare the above obtained equations for fl (or w), Eq. (1741) . with 
the equation of motion of 0. 

4.1.2 Comparison with the Goldstone picture 

The most general equation of motion for 0, Eq. (125!) . can be rewritten in the following 
form involving a total derivative 

V M |3A 5 V^0 + 3«V^ (d0) 2 + /3V" (0 ;(S7 ) 2 + 2/3V" (n# ; t) } = -^T. 

Hence, it is easy to see that, in the case of a spherically symmetric configuration 
<f)(R), this equation can always be integrated once to yield 

= i- [ R dRp(R] R\ 



R 3 R 3 Jo 

where we used that, according to Eq. ( 1671) . T/Mp = —p in the decoupling limit, and 
where C2 is an integration constant. Let us now see how this equation compares 
with equation (|69|) obtained for fl in the previous subsection. As we have seen, 
the relation between the unitary gauge (127|) and the gauge (1301) is encoded into 
the function p(R) though the coordinate change (I3TT) . On the other hand, it is also 
contained in the Goldstone field ir A , and in our case, in a single scalar field <fi defined 
as in ATMS} by 

n AB d B ^> = X A (x») - (79) 
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where as follows from (131!) X A and x M are given by = {t, x, y, z} and = 

{t, X,Y, Z}, and X A (x fM ) is given by the coordinate change (13"T1) . It is easy to see 
that equation (179)) has the solution (fi(R) 

<P' = d R (f> = R(e- l!i ^-l), (80) 



and that in this particular case no vector field A B is needed. Note that the relation 
between fi and is in fact "non perturbative" because of the presence of the ex- 
ponential in the above relation. The limiting conditions (134)) translate via relation 
(J80D into 

W(J2) = 0. (81) 

ix — >U 

This is precisely the limiting condition one would impose on a scalar field in a 
spherically symmetric configuration in order to have everywhere well defined second 
derivatives. For small fi, the above expression (1501) can be expanded as 

If we go to the canonically normalized 0, defined in equation (j!9p . and use the 
rescaled metric variables as defined in ([55]) . and then take the decoupling limit (1231) . 
one is left with the simple relation between and ft given by 

/* = • (82) 

Substituting 0' in the equation (1751) we get then 

3~ 1 l„ „ T -,~/2 . 3 D ~~// 1 d2 ~/~// 



and we recover exactly equation (169)) (identifying C2 and C2). It is however interest- 
ing to note that the boundary condition (15T|) can be lost in the decoupling limit due 
to the "non perturbative" (exact) relation between and fi given by equation (1501) . 
For example, if one considers the field theory defined by action ( I24jl and looks for a 
spherically symmetric solution, one is led to impose that 4>'(0) vanishes. This, via 
relation ( 152]) imposes however a condition on p, which is not necessary considering 
the conditions given at the end of section 13.11 
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4.1.3 Expected range of applicability of the Decoupling Limit 

We would like here to discuss to what extent the DL is capturing the leading be- 
haviour of the solution of the full nonlinear system. As we will see, the DL is 
expected to be valid for distances in-between the Compton length of graviton m^ 1 , 
and a scale that is parametrically lower than the Vainshtein radius and can even 
reach the Schwarzschild radius. 

First, as should be clear from the discussion of section [3.2.31 the DL should give 
a good description of the solution at least in the range Ry < R <C mT x . Another 
way to check this is to look more carefully at the expressions (j2~T|) . Indeed, it is easy 
to check that the quadratic, cubic and quartic non derivative self interactions (l2~Tj) 
(with appropriate values of ^1,^2,^3,^4) do not give significant^ corrections to the 
DL solutions for Ry < R <C m , while the cubic derivative <fi self interaction is 
retained in the DL. Comparing the other interactions fl2Tj) to the kinetic terms in 
the regime where the linearized theory is expected to hold (i.e. for R ^> -Ry), one 
sees that those interactions have each their own "Vainshtein" radii (corresponding to 
distances where those interactions generate 0(1) correction to the linearized theory) 
which are strictly smaller than Ry, as first noted in Ref. [IB]. Those radii are sent 
to zero in the DL, while Ry is kept unchanged. 

Let us now discuss what is going on at distances smaller than Ry. Using again the 
expressions ( I2T1) . it is easy to see that the quadratic, cubic and quartic non derivative 
self interactions stay negligeable at least up to R$, while one can estimate when the 
other interactions become of the same order as the cubic interaction retained in the 
DL. Indeed, the Vainshtein scaling (I61H63P translates into the scaling h ~ Rs/R, 
y/Rs/R (via Eq. flED])), and A ~ 0. From which one easily obtains the 
scaling of the canonically normalized fields h, <p, and A. Inserting this into Eq. (|2"T|) 
we see that a generic interaction (|2"T|) is much smaller than the cubic all the way 
down to the Schwarzschild radius Rs- Hence, as discussed in Ref. [16], assuming the 
Vainshtein scaling, one expects the DL to be correctly describing the solution of the 
full nonlinear theory below rrT 1 and down to Rs- The same question can be asked 
for the other scaling introduced above in Eq. (I6"l"|) . For this scaling, and below the 
Vainshtein radius, one still have h ~ Rs/R, while we have dd<p ~ ^ ~ m 2 RyR~ 2 . 
To estimate the range of validity of the leading behaviour of the DL solution, one 
has to pay attention to the fact that the leading term in Eq. (|64|) is a zero mode 
of the kinetic operator associated with the cubic term in the action of the DL. This 
kinetic operator, evaluated on this leading behaviour, is a sum of three terms which 
add to zero (see e.g. Eq. (173]) ). Hence, one expects to see the solution leave the 
found DL behaviour whenever the left over interactions (|2"Tj) generates a term in the 
equation of motion larger than, or of the order of, any of those three terms. For a 
generic interaction fl2~T]) this happens below the radius R^feM defined by 

Rktteto = Rv(Rsm)s + sW^ . (83) 
Again we see that this scale is always parametrically smaller than the Vainshtein 
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Radius, whenever k% + 2/c 3 — 6 is stricly positive. When this later condition is 
not fulfilled, it is also easy to see that the neglected interactions are then always 
subdominant with respect to the cubic self interaction retained by the DL for R < 
Rv- Hence, for the scaling (f6"4"j) . we expect as well the DL to give the leading 
behaviour of the solution of the full nonlinear system (if it exists) for distances 
below the Vainshtein radius and above the largest of the scales Rfaksfa which is 
simply R v x a. 

4.2 Solving the decoupling limit equation of motion 

4.2.1 Behaviour and boundary conditions at infinity 

As we have just demonstrated, in the decoupling limit, one is left with the single non- 
linear differential equation (175"]) to be solved with appropriate boundary conditions. 
Note, that in the simple case of BD and AGS potentials (respectively potentials (jl]) 
and (J5J), equation (175]) (taking then the simpler form (1771) ) can be put in the well 
studied form (with obvious notation), encompassing in particular the six Painleve 
transcendents 

w = F(w, £)w 2 + G(w, £)w + H(w, £)■ 

However, one can easily check that the functions F, G, and H above are such that 
the movable singularities of this equation are not only polar and even in this case 
solutions are not known analytically, hence one has to use a numerical integration 
in order to solve the equation of motion (1751) . 

This equation being of second order, we have to specify two initial conditions, 
w(C,i) and w(^i) at some point £j to start numerical integration from there. In our 
case, however, it is natural to start the integration from infinity, since we know the 
asymptotic behaviour of w there. Indeed, at infinity, i.e. large £ (or large R) regime, 
we are looking for a solution of equation (175*]) which has the asymptotic behaviour 
given by dropping all the nonlinearities. Hence it should behave at large £ as 

~ ™°o(£), (84) 

where u>oo(0 is defined by 

«U0 = gH- (85) 

This behaviour can then be used to integrate inward the equation of motion starting 
from the initial conditions 

= Woc(&), (86) 

Taking a sufficiently large & ^> 1, we may hope that a solution we obtain by 
numerical integration will be close to the true solution (if it exists) of Eq. (1751) with 
the asymptotic behaviour tu(£) — ► Woo(£,) as £ — > oo. 
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Things are slightly simpler in the case of equation (ITT)) corresponding to the BD 
and AGS potentials (0J and (J5]). Indeed, in this case, the change of variable 

Soo = r 3 , (87) 
= [w(0f 2 , (88) 

puts the differential equation (ITTj) in the form 

W 1/3 1 

<"^ + -» = 0. (89) 

In this form, the first term in the left hand side above stands for all the term 
quadratic in w in equation (!75|) (i.e the terms in Q), the second term represent the 
term linear in w and the last term encodes the source. The boundary conditions 
translate into the Cauchy initial values at 



^oo(O) = 

WL(0) = 0, (90) 

while the asymptotic behaviour at Soo — > is obtained now by dropping the second 
derivative of W m in the above equation (l89j) and reads 

\ 3/2 



(3 So0 ) 



In this simpler case, one sees clearly a crucial property of our problem, that holds 
also in the most general case (1751) . namely the Cauchy problem (I89H90I) (or (!75|) 
together with the asymptotic behaviour (184"]) ) is a singular Cauchy problem, since 
the differential operator of equation (1891) is obviously singular in (Woo, ^00) = (0, 0). 
As a result, one can not use standard theorems on Cauchy problem to conclude 
anything on the existence of a solution. Before turning to discuss solutions of our 
equation (J75J), we first discuss in the following subsection a simple linear differential 
equation with a similar singular Cauchy problem to illustrate some crucial properties 
also found in our solving of equation (1751) . 



4.2.2 A simple singular Cauchy problem 

Consider the following second-order differential equations: 

y"(x)+y(x) = - J (91) 

x 

-y"(x)+y(x) = -. (92) 

x 

These equations are singular at infinity similarly to Eq. (1751) . The term y(x) and 
1/x correspond respectively to the term linear in w and the source term l/£ 3 , while 
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the second derivative of y can be related to the nonlinear differential operator Q. 
Aside from the fact the differential operator of equation (l9"Tj) and (I9"2"j) is linear, those 
equation have the same structure as equation (j75J). The solution of (I9"T1) is given by 

yi(x) = G\ cos(x) + C 2 sin(x) + Ci(x) sin(x) — Si(x) cos(x), (93) 

and the solution of (I9"2"j) is 

y 2 ( x ) = C ie x + C 2 e- X - l - (e x Ei(-x) - e~ x Ei(x)) . (94) 

In Eq. ( 1931) and (I94"|) . C\ and C 2 are arbitrary constants, and Si(x), Ci(x) and Ei(a;) 
are sine integral, cosine integral and exponential integral functions correspondingly, 
given by 

«./ n Tsint, . . pcost, , /' 00 exp(-t), 

Si(x) = / dt, Ci(x) = - dt, Ei(x) = - — — — ^dt. 

JO t Jx t J_ x t 

Now let us find solution(s) of (19"T|) and (I9"2"j) such that, 

2/i 2(2?) - * — , when x — > 00. (95) 

This correspond, mutatis mutandis, to the asymptotic at infinity, Eq. (!85j) . There 
is a unique solution of (|9"Ii) with the asymptotic f l95l ). 



7T 

= — cos(x) + Ci(x) sin(x) — Si(x) cos(x), (96) 

while for the equation (19"2"1) the asymptotic behaviour (|95p does not fix uniquely the 
solution, 

y 2 (x) = C 2 e~ x - X - (e x Ei{-x) - e~ x Ei(x)) . (97) 

I.e. one has the freedom to choose at will the integration constant C 2 . Solution (196|) 
has the following asymptotic expansion 

, x 1 2 24 720 n ( 1 

yi(z) = 3 + — - — + — 

tXj JU JU JU \ JU 

while looking for a power serie expansion of solution (1971) . we find 

1 2 24 720 ^ / 1 

y 2 (x) = - + — + — + — + — 

X X 6 x° x' \X J 

Those series expansions are in fact divergent, and hence only give at best asymptotic 
expansions of the solution, however we notice that the last one misses the existence 
of the homogeneous mode C 2 e~ x . One can nonetheless easily see that the solutions 
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(196]) and (l9Tj) are in fact finite and they give the correct asymptotic (195]) at the 
infinity. 

The toy examples (1911) and (I92p show that for singular equations, the asymptotic 
behaviour at infinity can be enough to fix uniquely the solution (and hence the 
singular Cauchy problem is well-posed), this is the case for example of Eq. (19~T]) ; 
while in other cases, such as that of Eq. (192]) . there exist a family of solutions with 
the same asymptotic behaviour, which can be missed by looking for the solution 
by a power serie expansion. The properties summarized in this last paragraph are 
exactly recovered for our equation (!75l) . as we now see, beginning by discussing the 
simplest case of BD and AGS potentials. 



4.2.3 Solutions for the BD and AGS potentials outside the source 

In the case of the BD and AGS potentials (of Eq. (J5J) and §5§ respectively), the 
nonlinear differential equation to solve takes the simple form (!77|) . In those cases, it 
is easy to find a power serie expansion of the solution around £ = +00 in the form 

00 

»(« = E^k- (*) 

n=0 ^ 

The first non trivial six w n coefficients, are given by 

2 4 1024 712960 104910848 225030664192 
~ W + S W + W* +S 24H [E+ 243 e 3 +S 2187 e 28 + "' ( } 



Note that the first two terms of this series match those given in equation (1571) . and 
in fact this expansion is just the continuation of the one given in this equation. 
One can check numerically that this power series is in fact divergent. However we 
found it useful for the purpose of numerical integration to use this expansion as an 
asymptotic expansion of the solutions. Indeed it provides a better approximation of 
the solution than keeping only the leading order (186]) . provided that we truncate the 
expansion at the appropriate order 0- In both cases of BD and AGS potentials, we 
were able to integrate numerically inwards equation (1771) and obtain non singular 
solutions all the way to small radii. Our results are shown in figures [1] and [3) 



17 I.e., as usual with asymptotic expansions, for any given £j from which we want to integrate 
inwards, there is an order n that minimizes the n-th term w n /^ 3+5n of the series (|98[) . and we 
found it helpful to use the formal series truncated at this order n«), i.e. to replace the initial 
conditions ([86]) by 



= E 



Wk 



£3+5fc ' 
k=Q S 

(„(«,) 
E ^ I • (100) 
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Figure 1: Solution for w in the case of the BD potential. The numerical solution 
is shown by solid thick (blue) line. For distances much larger than the Vainshtein 
radius, £ ^> 1, the solution is well approximated by the asymptotic fl85|) . shown 
by dashed thin (black) line. Close to the source, ^ < 1, the numerical solution 
approaches the Q-scaling asymptotic, Eq. (11031) . shown by dotted thin (black) line. 



These numerical solutions exhibit different behaviours at small £ which can be 
understood as we now explain. First, if, following Vainshtein, we assume that we 
can drop the term linear in w in equation ( 1771) . the leading behaviour should be 
given by a solution of the equation 

2QW = ^, 

namely the same equation as equation fl59l) . An assumed power law behaviour leads 
then to the Vainshtein scaling at small £, as (see Eq. fl63l) ) 

1 



W ~ Wy oc 



(101) 



Whether such a scaling leads to a real or an imaginary solution depends on the 

potential. For the AGS potential it leads to a real leading behaviour w ~ wy = yj^- 

An expansion of the solution around this Vainshtein scaling can then be obtained 
in the following form 



»(0 = ^ 




B C 



3x/s 3 (—5 + v^) d2 e _ 2 +^I , 6£ 2 



8^ (-4 + 



3^-^ + ^ + - (102) 
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This expansion can be seen to match correctly the (only - see thereafter) solution of 
the AGS case that has the small £ Vainshtein asymptotics, provided the constant B 
is fixed to some specific value. It can be understood as follows. The leading term is 
of course the Vainshtein scaling. The next term oc ^i(^ J+3 ^) is a zero mode of the 
linearized Q around the leading behaviour. The rest of the series is then obtained, 
as usual, order by order in £. 

For the BD potential, however, the Vainshtein scaling at small £ leads to an imag- 
inary solution, which is not acceptable. In the BD case, rather than the Vainshtein 
scaling, we found that there exists a real non singular numerical solution all the way 
to small £ which interpolates between the large distance behaviour w ~ l/£ 3 for 
£ 3> 1 and the small distance behaviour w ~ l/£ 2 for £ <C 1 (see Figured]). Notice 
that this solution is perfectly regular around the Vainshtein radius £y = 1. The 
found new scaling at small £, namely 

w ~ i/e, (103) 

and named in the following Q -scaling, is in fact obtained from a zero mode of 
the operator Q. Indeed one can check that the equation Q(w) = has the exact 
solution^ 

K \* 



w q(0 = \^ + K 



C 3 



2 /3 1 /3 

where K and K\ are two arbitrary constants and A = K ,B = {2Ki)/(3K ) 
are the first two coefficients of the expansion of the zero mode around £ = 0. A 
solution of equation (1771) can then be found in the form of the double expansion 

. oo n 

^) = ^+EE w »^>o fc (104) 

^ n=l fc=0 
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This solution can be easily found using the new function and variables 

6 = s 1/3 

which translates the differential equation (|77|) into the following nonlinear equation on Wo (So) 

jfW Q {t) W 1/3 1 

dBg 4^ 3 + 6S < /3 ' 
where Q is transformed into the first term on the left hand side. 
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where the coefficients w n ^ are functions of A 0} B . For the BD and AGS potentials, 
the first coefficients read: 



"(0 =f ^ 



A 3A B - s lnf 



3A 



1 + s 6A Bq - 54A 2 £? 2 - (2 - s 36A B ) ln£ - 6 In 2 g 4 5 

+ 216A 3 , * + ^ j ' 

(105) 

The values of the parameters A Q , B can be fixed matching numerically the expansion 
at small £ with the large £ asymptotic behaviour. In the case of the BD potential, 
we find A ~ 0.645 and B ~ 0.208 and we find that the expansion fits very well 
the numerical solution up to £ ~ 1 (the Vainshtein radius), as can be seen in Fig. [21 
Expansions for w(£) and v (£) can be found in the same way. They read 
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Figure 2: Plot of the numerical solution (solid black curve), and the series expansion 
(11051) given up to order £ (dotted red curve) and up to order £ 13 (dashed blue curve). 
The expansion up to £ 13 approximates the numerical solution with a precision better 
that 99% in the range ^ £ ^ 0.95. 
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«(0 = \-Y~& {3AoBo ~ sln °~^ + ^ (106) 

«(0 = -~ + (A) + ~A,ln£) 

+^ (s + 9A i? - 5 3 InO + + O (£ 5 ) . (107) 

In dimensionful units, and keeping only the dominant terms, the above equations 
(11061) and (11071) read (the corresponding expression for \x has been given in equation 

AMD) 

A = ^ + 0(1) (108) 
v = -^-^ l -A,{R s MflHnR + 0{l) (109) 

Note that the leading terms at small R, A = Rs/R and v = —Rs/R, correspond 
to the linearization of the Schwarzschild solution in General Relativity, and hence, 
even with our new scaling, one recovers (in the decoupling limit) the Schwarzschild 
solution at distances smaller than the Vainshtein radius, in agreement with the 
Vainshtein original idea that nonlinearities can cure the vDVZ discontinuity. Notice 
however that this happens here for a potential, the BD potential, for which the 
Vainshtein mechanism was believed not to work because of the imaginary nature of 
the Vainshtein scaling at small £ [12J. It is also interesting to notice that the first 
correction to the Newtonian potential arising from equation (11091) has a form similar 
to the one required by MOND pp. 

Our numerical investigations of the solution in the BD case lead to the con- 
clusion that there is in this case a unique solution (shown in figure [T]) with the 
right asymptotics (1851) . This was confirmed by an analytic proof of this uniqueness 
(around £ = +oo) provided by J. Ecalle [28]. To conclude on the BD potential, let 
us further mention that in this case, the numerical integration appears to be stable 
against small perturbations. Namely, a slight change in the initial conditions at large 
£i does not result in growing divergences but only in small oscillations around the 
found unique solution. This can be understood analytically by perturbing around 
the asymptotic solution (see appendix ICl) . 

Let us now turn to discuss the AGS case (potential (jSJ), or s = +1). In this case, 
the expansion (|99~I) still holds, but in analogy with the example (|92|) we found in this 
case infinitely many solutions with the same power serie expansion (199]) at infinity. 
In this case, indeed, we found that with the same asymptotics at infinity, different 
behaviours, including the original Vainshtein scaling, were possible at small £. In 
fact, the most general behaviour at small £, is still of the form of the Q-scaling (I103[) 
with the asymptotic expansion (U04H105p . but depending on the chosen solution at 
infinity, this behaviour can be picked up at arbitrarily small £ while there is an 
intermediate region at £ smaller than one where the solution follows the Vainshtein 
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scaling fllOip (see Fig. [3]). The Vainshtein scaling then appears as the limiting 
case of the family of solution having the right asymptotics as infinity (1851) . and 
so can be obtained in the AGS case. From a practical point of view, numerical 
integration in the AGS case is quite difficult, due to numerical instabilities which 
appear for distances larger that the Vainshtein radius £ = 1. Indeed, any small 
departure from the solution sources exponentially growing modes and Runge-Kutta 
type of integrations reach singularities quickly. In this context, we found very useful, 
while integrating from large £j towards £ = 0, to use the asymptotic expansion flUSl) . 
truncated at the order appropriate to £j. We were also able to confirm our numerical 
results by solving equation (175|) using a relaxation method, which is more robust 
under numerical instabilities. Again, our numerical investigations were confirmed by 
mathematical proofs of the existence (around £ = +oo) of infinitely many solutions 
with asymptotics (1851) [28, 29J. In fact, one can also show that any two solutions 
of (IToT) having the same asymptotics (185]) differ (at dominant order) at large £ by a 
quantity given by a constant times £ 3 / 2 exp (—k 3/5 £ 5 ^ 2 ) where k is an integer [28J. 
Finally, one can also analyse the problem by perturbation theory around the large £ 
behaviour (I55"j) . This is done in appendixOand we find that there are both a growing 
and a decaying mode. The growing mode is responsible for the numerical instability 
pointed out above and should be discarded in the true solution; the subdominant 
decaying mode, similar to the decreasing exponential found for equation (|92j) . can 
be freely specified while keeping the asymptotic behaviour (!85|) . 



4.2.4 Solution for the BD and AGS potentials inside the source 

Let's now include a source. For simplicity, we consider a star of constant density, 
p a , and radius £ (i.e. of radius R Q in physical units). We have 

3 

Pa — 73-, 

leading to the equation for w inside the source 

~, s 3 1 
2Q(w) + -w = w . 

Here too, two different types of solutions can be guessed. The first kind are solutions 
of the form 

A °° 

«(0 = i£ + ( n °) 

^ 71=1 

corresponding again to a l/£ 2 leading behaviour for small £. The second kind are 
solutions of the form 

mo = $>^ 2n = a + s 6 -^f-^ e + ... (in) 
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Figure 3: Plot of the numerical solutions for w in the case of the AGS potential, 
starting from £j = 2.5. For £ ^> 1 all the solutions approach the asymptotic fl85l) . 
shown by dashed thin (black) line. At small distances, ^ < 1 the solutions pick 
up different asymptotic regimes. The solid thick (blue) line corresponds to the 
Vainshtein scaling w ~ l/VC- The dash-dotted thick (green) line corresponds to 
a solution with the Q-scaling (11031) . The dashed thick (red) line corresponds to 
a solution which first follows the Vainshtein scaling and then finally picks up the 
Q-scaling (TTOBl . 



where the leading term is now w ~ A$ (and s = — 1 in the BD case, and s = +1 in 
the AGS one). 

Our numerical integration indicates that, the BD solution, which is fixed by the 
asymptotic condition at infinity (|86|) . follows the w ~ l/£ 2 behaviour inside the 
star, as shown in Fig. HI In the AGS case, the situation is more subtle: generically 
(i.e. without fine tuning), the non singular solutions also adopt such a l/£ 2 leading 
behaviour; an example of such a solution is presented in Fig. (thick dash-dotted 
red curve). However, there is a case for which the other scaling w ~ A is possible. 
This case precisely corresponds to the Vainshtein solution outside of the source, as 
shown in Fig. [5] (thick blue curve). 

Note that for the Vainshtein and Q-scaling at small R, and in the absence of a 
source, the Jacobian ( 1321) is going asymptotically to zero. This is not true for the 
scaling f ill 11) . Moreover, only this scaling leads to a finite Ricci scalar of the physical 
metric at R = 0. Hence, this might provide a way to select a physical solution in 
the full nonlinear theory. Discussion of these interesting issues is however left for a 
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Figure 4: Plot of the numerical solution for w in the case of the BD potential in 
the presence of a source of radius £ = 0.01 (depicted above by a light blue area), 
starting from £j = 6. The blue thick line corresponds to the Q-scaling of Eq. ( 11 101) . 
which is picked up by the asymptotic behaviour at infinity. Note that the solution 
is almost not affected by the presence of the source. 



future work [30] , since it involves keeping track of nonlinear terms in the equations 
beyond those appearing in the DL. 

4.2.5 Solution in the most general case 

In the most general case of equation (I75p with arbitrary a and (3, a solution with 
a Vainshtein scaling at small £ can easily be found. It is obtained by a power law 
ansatz and dropping the term linear in w in Eq. (l75|) and reads 



This solution is real iff 
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3(25a + 13(3) _ 



1/2 



-1/2 



as £ — > 0. 



25a + 13(3 > 0. 

In analogy with the discussion of the previous subsection, one can look for another 
scaling where the dominant term (at small £) is given by a zero mode of the operator 
Q appearing in equation ( 1751) . Zero modes in the form of power law, 



WQ = A? 



(112) 
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Figure 5: Plot of the numerical solutions for w in the case of the AGS potential 
in the presence of a source of radius £ = 0.01 (depicted above by a light blue 
area), starting from ^ = 3.5. The solid thick (blue) line corresponds to the solution 
having the behaviour Eq. (11 111) inside the source. The thick dash-dotted (red) line 
illustrates the Q-scaling, Eq. (lllOp . inside the source. 



can be found, where p is given by 

P± = r4 > for a + /?^0, 113 

a + p 

where the plus sign in the right hand side corresponds to the power p + , and the 
minus sign to the power p_. Note that for a + (3 = 0, e.g. the BD and AGS 
potentials studied in subsection 14.2.31 and 14.2.41 the Q-scaling p = —2 can be found 
from (I113p taking the limit a —>■ —f3. One can check also that, in order to describe 
the dominant behaviour at small £, any particular root, p + or p_, should satisfy an 
additional criterion, namely, 

1 

y 2' 

otherwise with this scaling the nonlinear terms in Eq. (175]) are not dominant. 

The various possibilities for the different power law scaling as £ — >• are summa- 
rized in Fig. El in the (a, 0) plane. Roughly speaking Eq. (j75|) possesses four different 
regimes at £ — > with: 

• only Vainshtein scaling; 
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Vainshtein scailim 



Q-scaling 



both scalings 



Figure 6: The diagram of possible scalings for solutions of Eq. ( 1751) at £ — > in the 
plane of parameters a and f3. The whole parameter space (excluding the origin) can 
be reduced to the circle of unit radius, with the help of an appropriate rescaling (see 
the text). 



• only Q-scaling; 

• both Vainshtein and Q-scalings; 

• no power law scaling. 

There are three boundary lines, where interesting transitions from one type to an- 
other can happen: 

•0 = 0; 

• a = -13/3/25; 

• a = -(3/2. 

Note, that with the help of the rescaling £ — > and w — > C~ 3 w, with C = 
(e* 2 + /3 2 ) 1/10 , Eq. m 

can be brought to the same form where the coefficients a 
and (3 now verify a 2 + 1 = 1, and hence one can restrict the study of the different 
cases in the (a, (3) plane of Fig. [6] to those where a and (3 lie on the circle of 
unit radius. This circle is mapped to the horizontal segment in Fig. [7J which shows 
together with the possible power scaling (the values of the Vainshtein and Q-scaling 
power scalings are plotted on the vertical axis), results of the numerical integrations 
of equation (175p indicating which of the scalings are in fact realized in the solutions 
having the right large £ asymptotics given by (1851) . Depending on the sign of 
(3, this asymptotics does or does not fix a unique behaviour at small £. Indeed, for 
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Figure 7: The values of the possible power law scaling exponents of the solutions of 
Eq. (|74|) at small £ as a function of parameters a and f3. The circle in the Fig. [H] 
is mapped to the horizontal segment. The Vainshtein scaling is shown by dashed 
horizontal (red) line, while the Q-scalings are shown by solid (green) lines. The 
scalings which are found to smoothly continue to the asymptotic (1851) at large £ 
are shown by thick lines. The scalings at £ — > which do not give the required 
asymptotic at large £ are shown by thin lines. 

(3 > small variations of initial conditions at large £ does not affect the solution at 
small £, similar to the case we have already studied, the BD potential. This signals 
about the uniqueness of the solution for Eq. (|75p when the asymptotic behaviour is 
fixed by flSBT) . On the contrary, for (3 < 0, the asymptotic behaviour flHB^) does not 
fix uniquely the solution, and we find infinitely many different solutions at small £ 
with the same behaviour at large £. 

Starting from the point (a, (3) = (—1, 0), and moving counter-counter-clock- wise 
along the unit circle of figure [HI i.e. from the left and along the horizontal segment 
of figure [3, let us now describe with more details the successive cases of interest 
together with the results of our numerical integrations. 

• (3 = 0, a < 0. 

This point is somewhat special: no Vainshtein scaling is possible and the two 
Q-scalings are equal and coincide with the asymptotic solution at infinity, 

In fact, in this case the asymptotic form Eq. (|85|) satisfies Eq. (!75|) exactly for 
any £ > 0, 
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This can be easily understood recalling that whenever (3 vanishes, the covariant 
form of Q^ a '^ can be read off in equation (1251) to contain only □ acting on <p, 
and hence, in vacuum, a solution of □</> = is a solution everywhere^!. 

P > 0, -oo < a < -(3/2. 

The solution at £ — >• 0, with the correct large distance asymptotics (1851) . is 
always given by the Q-scaling with p_, although in some parts of this region 
of parameters we could expect the second Q-scaling and/or Vainshtein scaling. 
This region contains the BD potential, which we have already studied in detail. 

P>0, -P/2 <a< +oo. 

The solution is always singular in this region, the integration from large £ 
always breaks down at some finite distance £ > 0. One could have expected 
the existence of a regular solution with the Vainshtein asymptotic scaling but 
it appears that this scaling does not give the required asymptotic behaviour 
at infinity. 

P = 0, a > 0. 

This case is similar to the case, P = 0, a < 0. An exact solution for £ > is 
given by an analytic expression, 

p <0, ~13p/25 <a< +oo. 

For p < the asymptotic condition at infinity (|85|) does not fix uniquely the 
solution, as in the case of AGS potential, which is included in this region 
of parameters. Thus tuning the initial conditions at large £ we are able to 
obtain at small £ either the Q-scaling with pJ^ or the Vainshtein scaling or a 
singular solution. It is remarkable that the other Q-scaling cannot be obtained 
by the tuning of the initial conditions even though the scaling with p_ would 
be a dominant scaling for — P < a < +oo. Here again it is found, as in the 
AGS case, that the Vainshtein scaling appears as a limiting regime of a family 
of solutions which all have the Q-scaling as their innermost scaling, but can 
behave at some intermediate range of £ following the Vainshtein scaling (see 
figure [3]) . 

P <0, -P/2 <a< -Up/25. 

The only possible scaling here, as confirmed by numerical integration^!], is a 
family of the Q-scaling with p + and different A in fll 12j) . Of course, changes 
in the initial conditions also lead to singular solutions as well. 



19 We leave for a future publication the discussion of the inclusion of a source. 
20 Note that this is a whole family of solutions, since by tuning the initial conditions at large £ 
we can get different factor A in (|112p at small £. 
21 the border case with p + = p- is included. 
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• (3 < 0, -oo < a < -(5/2. 

There is no scaling possible for small £, and indeed, the numerical integration 
gives only singular solutions in this region. 

When it is possible to find a scaling for small £ (either of Vainshtein or of Q type), 
a full solution can be obtained as a power law expansion in £, similar to the ones 
obtained in the specific cases of AGS and BD potentials (Eq. (11021) and (1 1 5 [) ) . Note 
that for both Vainshtein and Q-scalings, it is always possible to add to the leading 
behaviour a zero mode of the linearized Q around the leading term, introducing a 
free constant (denoted B in Eq. (11021) and f 1 1 5 p ) which should be fixed by matching 
the small £ expansion with the asymptotic behaviour (1851) . 

To conclude, we have found that for the range of parameters (3 < 0, a < —13/2 
and (3 > 0, a > —(3/2, regular solution(s) with the right large distance asymptotics 
fl85|) exist. For positive (3, and fixed a only one such solution exists and it has the 
Q-scaling at small £. For negative (3, and fixed a, a whole family of such solutions 
exist, which have the Q-scaling at small £, while the Vainshtein scaling appears as a 
limiting case of this family of solutions. As such, and in contrast to the Q-scaling, it 
never appears as the only possible scaling realized at small £. In the other range of 
parameters only singular solution(s) can be found if the asymptotic behaviour (|85|) 
is fixed. 

5 Discussion and conclusions 

In this paper we have investigated static spherically symmetric solutions of nonlinear 
massive gravities in the so-called decoupling limit (DL). After having first identified 
how to obtain this limit with the ansatze we used and which are appropriate for 
studying static spherical symmetry, we have solved the system of equations left over 
in the DL and obtained non singular solutions featuring a Vainshtein-like recovery 
of solutions of General Relativity (GR). This first shows that the singularities found 
to arise solving the full nonlinear system of equations [12] are not present in the DL, 
despite the fact those singularities are commonly thought to be due to a negative 
energy mode also seen in this limit and associated with the kind of instability first 
discussed by Boulware and Deser [11]. To us, this is not necessarily a surprise, 
since there is in fact no clear clash between having (static), non singular, spherically 
symmetric solutions and a ghost (coupled to positive energy modes) in a model. 
In fact, appendix [D] gives a toy example with such a property. One should of 
course worry about the possible instabilities (in time) of the static solutions, but 
this is another story and has a priori nothing to do with the singularities discussed 
in [12]. Moreover, we also found that the scaling at distances smaller than the 
Vainshtein radius, first conjectured by Vainshtein JTU], was only a limiting case in 
an infinite family of non singular solutions each showing a Vainshtein recovery of 
GR solutions below the Vainshtein radius but a different (but common) scaling at 
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small distances. This new scaling was shown to be associated with zero modes of the 
nonlinearities left over in the DL. Interestingly, this family of solution all have the 
same asymptotic power law expansions but differ at small distances which means in 
particular that the asymptotic expansion is not enough to fix uniquely the solution 
at small distance. We also noticed, including sources, that the Vainshtein scaling 
gives a better behaviour at the origin, but it is hard to conclude on this matter 
without studying the nonlinearities not included in the DL and we plan to come 
back on this issue in a future work [5U]. Last, we have also shown that potentials 
that were thought not to admit a Vainshtein-like recovery of GR solutions, because 
the conjectured Vainshtein scaling at small distances would have lead to imaginary 
solutions, could in fact accommodate such a recovery via a different scaling at small 
distances. This scaling is the same as the one mentioned above associated with zero 
modes of the nonlinearities appearing in the DL. 

An intriguing question is whether the properties found in the DL, and summa- 
rized above, also hold in the full nonlinear system. This requires in particular solving 
numerically the full nonlinear system, which is not an easy task given the numerical 
instabilities we already noticed in the DL. We will carry a thorough analysis of the 
full nonlinear system, stressing in particular the role of the neglected interactions, in 
a future publication [SU]. At this stage, let us just say that the findings of this work 
opens the possibility that important properties of the full nonlinear system could 
have been overlooked. On the other hand, if indeed this is not the case, it shows 
that the DL is not capturing all the interesting physics of the spherically symmetric 
solutions. It is also important to answer this question for the sake of a better under- 
standing of other models, such as the DGP model where the role of nonlinearities 
in the equivalent DL and beyond are not fully understood. E.g. there are no known 
non singular exact solution in the bulk and on the brane describing the equivalent of 
a Schwarzschild Black Hole in DGP theory and results on the Vainshtein mechanism 
in DGP gravity are all obtained using some approximation scheme [T4"l [TBI EI]- 
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A Equations of motion of X 



A 



We first derive the form of the variation Sf^x) of /^(x), defined as in equation 
(TTTj) . when X A varies as X A {x) — > X A (x) + 5X A (x). We have, replacing X A (x) by 
X A (x) + 5X A (x) in the definition of f fll/ (x), that 5f fll> (x) is given by 

SMx) = d,5X A d u X B f AB (X) 
+d^X A dJX B f AB (X) 

+d ll X A d l/ X B 5X c d c f AB (X). (114) 

Using then that 

8 G A = d A x°d a X c , 

where x a (X) stands for the inverse mapping of X (x), we find that 

d v X B f AB {X{x)) = dAx'Uix), 

and 

dcfAB(X) = d c x°d a X D d D f AB {X). 
Inserting those expression in (j!14p . we find that the latter reads 

Sfp, = d,5X A d A x°f au (x) 
+d u 5X B d B x°f^(x) 

+5X c d c x°d (7 X D d ll X A d u X B d D f AB (X). 

We then define 8x tl {x) as in (fl3l) . and rewrite df^ as 

5fnv = d^5x a f au + d u 5x a f^ + 5x a d a f fMU 

-8X c f EF d v X F (d B d c x°d,X B d a X E + dcx'd^X*) 
-8X c f EF d^X E \d A d c ^d v X A d a X F + d c x°d a d v X F ) (115) 

If one differentiates with respect to x M (respectively x v ) the quantity 

8 E = d c x°d a X E , 

we see that the last terms in the last two lines of the equation (11151) vanish, and hence 
we get that 8f^ u is given by the expression (fT2l) . that is to say the Lie derivative of 
fnv along the expression 8x^ considered as a vector field on the space-time manifold. 
In fact it is easy to see that 8x^ transforms as a vector field under a coordinate 
change x'^ = x /fl (x). Indeed, under such a coordinate change X A and 5X A are 
scalar quantities, and hence we get that 

8x'Hx') = 8x a (x(x'))^—, 

ax a 
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using the fact that x'^{X'(x')) which appears in the definition of Sx >fl (x') can also 
be rewritten as 

x'»{X'{x')) = x'»(x(X(x(x')))). 

We now use the fact that the interaction term between the two metrics, appearing 
in the action <Si nt is a scalar under reparametrization. Hence, under a coordinate 
change of the form x^ — > x M + one has that 

(^ v< " ,(/ ' 9) ) £!/ - + (^M v, ° ,(/ ' 9, ) £ ^ 

is a total derivative, so its integral over space-time vanishes (provided we consider 
asymptotically vanishing £ fields). Using then that 

and the definition of the energy momentum tensor (jUj), we obtain that under a 
reparametrization generated by the infinitesimal vector field one has 

J d 4 xV^~9i^ u T^ = - l -m 2 M 2 P J d'x ^J^vM{f,g)) C^. 

Applying this formula with £^ = 5x^ defined above, we get the expression ffl4|) . 



B Quadratic part of the equations of motion in 
the A, //, v gauge 



We find for Gi? and Ggg 

while the quantities Q tt , Qrr and Qb are depending on the choice of interaction 
term Si n t[f, g] in action ([3]). E.g. for the two interaction terms (Tj0) and (111), they 
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read respectively 
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As we now show, in full generality Q(fi) depends only on two parameters, a and /3, 
leading to the expression of Eq. fl54|) . We start from the general expression for the 
interaction term ([6]) 



S, 



1 



int 



m 2 M 2 P 



d 4 x V(g,f) 



1 



m 2 M 2 P 



d A x^g V{zH). 



The scalar V A (g~ 1 f) depending only on four invariants made out of the rank two 
tensor M{? = ' j uv (or using a matrix notation M = g _1 f), we choose here those 
invariants A* to be given by 



Aj = tr [(M - ll) 1 ] 



1,2,3,4. 



Thus, one can write 



The stress-energy tensor then reads 



M 2 m 2 



1 T/ A 



5A, 



^ 5^ <9A 4 



where ^ = (M - 1)* f, using again a matrix notation. We can restrict now the 
metrics to be of the spherically symmetric form (J28~l) . Since we are interested only 
in the terms quadratic in /i, we can further specify our study to the case for which 
A = v = 0, and expand the stress-energy tensor to the second order in /t, which 
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requires to expand the potential V to third order in //. This can be done easily 
noticing that Aj = 0(fi l ), leading to the expansion: 

V(Ai) = ax + a 2 A 1 + a 3 Aj + a 4 A 2 + a 5 Af + a 6 AiA a + a 7 A 3 + C(/i 4 ). 

With this expansion, the Bianchi identity (1351) reads 

= (116) 

with Q (a >P\fi) defined as in Eq. (I54p . and 

= — a2 — 4a3 — 4a4 , 
a = — — (a 2 + 4a 3 + 16a 4 - 48a 5 + 24a 7 j, 

^ = -— ( - 7a 2 - 68a 3 + 88a 4 - 80a 6 - 120a 7 ) . 

The last step is to impose the Pauli-Fierz form ([1]) for the mass term, which fixes 
the terms linear in \x in the right hand side of the first line of Eq. (11161) to be zero, 
i.e. f2 = 0. Identifying then the remaining quadratic term Q^ a ^\jx) in fl 1 1 6 [) with 
the quadratic term Q((j,q) in Eq. (I5T1) leads eventually to the general expression of 
Eq. (j54"j) for the quadratic term in \i in the Bianchi identity. Hence, as said in the 
main body of the text, it is here explicitly seen that the fact no term linear in \x 
appears in the Bianchi identity, a crucial property for what concerns the scalings 
appearing in the Vainshtein mechanism, is due to the peculiar tensorial structure of 
the Pauli-Fierz mass term. 



C Expansion around £ = +00 

Here we study with some details solutions of Eq. (1751) for general a and /?, expanding 
around the leading behaviour (1551) . It is convenient to factorize the asymptotic form 
of w given by (1551) and further use the following change of parametrization 

w (0 = ^[i-e 3 G(0] 

2/5 



The asymptotic behaviour, (155|) w ~ -u^ = 2/(3£ 3 ), translates into 



ec(o - o 

ecG'io - o (H7) 
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when ( — > oo. In terms of the new variables ( and G, Eq. (!75l) takes the form 



5(0 



9a + 5/5 



+ 



1 - 



4/5 

11 9a + 5/3 
20 4/? 



e 3 G(o 



+ ^(C), 

with the source S(() being given by 



15 a + (3 
~8 6 



G"(() 



— ' 33 ^W'(o 



G'(C) 



118) 



e 



It is then clear that the limiting behaviour (11171) for G allows to drop all the nonlinear 
terms in G the equation flllSp in the vicinity of infinity. This procedure leads to a 
Bessel equation of 0-th order with a source: 



G"{Q 



G'(C) 

c 



+ l^G(C) = S{Q. 



;ii9) 



For (3 > 0, the general solution of (11191) is given by 



7T 



7T 



G(C) = 2 MO J Y (t)s(t)tdt - -y (C) jf Mt) s(t)tdt + a j (c) + d 2 y (Q ■ 

The only possible choice of integration constants Di and D 2 consistent with the 
conditions (11171) . is that they both vanish, i.e. (Di, D 2 ) = (0, 0). In this case, G has 
the asymptotic behaviour G(() = 0((+ 16 / 5 ) and the asymptotic conditions (11 lTj) 
fixes the solution uniquely. This is similar to the example (I9T1) discussed in the text 
and also allows to understand qualitatively properties of the numerical integration of 
equation (1751) for positive (3. Indeed, in this case, if one integrate numerically inwards 
(from large £), numerical errors can source the homogeneous modes «/o(C))^o(C)- 
However, the later are not growing out of control in agreement with the fact the 
numerical solution is found to be stable under such small perturbations. 
For (3 < 0, the general solution of (11191) reads 

/C r°° 
I (t)S(t)tdt - Jo(C) ^ K (t) S(t)tdt + D 3 K (() + ZVo(C) • 



It this case, the conditions (11171) only leads to the vanishing of _D 4 , but the other 
integration constant D 3 remains free, since the homogeneous mode i^o(C) decays 
fast enough at infinity. Whatever the chosen D 3 , one finds a leading behaviour for 
G given by G(() = 0(C 16 ^ 5 )- This is analogous to the example (|9~2"|) discussed 
in the main body of the text. The existence of the growing mode /o(C) allows us 
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to understand better the numerical instabilities observed in the case of the AGS 
potential. It is indeed impossible not to source this growing mode numerically, and, 
after some point, this mode dominates the solution and the integration reaches a 
singularity. As a consequence, one has to tune very carefully the initial conditions 
of the numerical integration in order to avoid this explosion. 

D An example of a system with coupled normal 
and ghost fields. 

Let us consider the following higher-order derivative action, S, for the scalar field <fi 

S = J d A x | ~ (D0) 2 + T<f>y (120) 

The above system can be described in terms of two fields, one positive energy ("nor- 
mal") field, <f) c , and one negative energy (ghost) field, tp c . In terms of these new 
variable, S now reads 

S = Jd A x |^ C D0 C - -\(<t>c- + ^(<t>c + V>c)} > (121) 

where it is explicitly seen that the normal and ghost degrees of freedom are coupled 
via a potential term. From (11201) we find the equation of motion for <fi, 

□ (□</>) =T. (122) 

In the case of the point-like source, the full solution of (11221) is, 

/ At R AnR? Bi „ 

4 > (R) = ^ + ^—-^ + B . (123) 

The first term in ffT23l "kills" the delta-function in the r.h.s. of ffT22l . Note, that the 
solution is smooth everywhere for R > 0. Moreover, the presence of the growing with 
R terms in (11231) is not dangerous, since these terms do not lead to a "bad" physical 
behaviour at R — > oo. Indeed, the energy density, E(R), for the configuration of (j> 
given by ffT23l) is [31] 

^)4(^) 2 +#(t+^+§)- 

Thus, E(R) -> (1/2)A 2 as R -> oo. 
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